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ABSTRACT 

Image computation is a fundamental tool for performance assessment of astronomical 
instrumentation, usually implemented by Fourier transform techniques. We review the 
numerical implementation, evaluating a direct implementation of the discrete Fourier 
transform (DFT) algorithm, compared with fast Fourier transform (FFT) tools. Sim- 
ulations show that the precision is quite comparable, but in the case investigated the 
computing performance is considerably higher for DFT than FFT. The application to 
image simulation for the mission Gaia and for Extremely Large Telescopes is discussed. 

Key words: telescopes - methods: numerical - instrumentation: miscellaneous. 



1 INTRODUCTION 

The imaging performance of astronomical instruments, and 
optical instrumentation in general, is related to the Point 
Spread Function (PSF), i.e. the intensity distribution from 
a point-like object at infinity, which is in most cases a sat- 
isfactory representation of a star. The monochromatic PSF 
is computed as square modulus of the field amplitude on 
the focal plane; the amplitude function, in turn, is built 
through the diffraction integral, derived from the wave de- 
scription of electromagnetic radiation. An historical review 
of the diffraction theory is out of the scope of this paper, de- 
voted to performance comparison of two Fourier transform 
(FT) algorithms computing the diffraction integral. 

The PSF includes, by a suitable wavefront description, 
details of the instrumental response. The polychromatic PSF 
of an unresolved star is produced by superposition, weighed 
by spectral distribution, of the monochromatic PSFs; thus, 
objects with different spectral characteristics have different 
polychromatic PSF. The image of a resolved source can be 
described as a superposition of PSFs, weighed by spectrum 
and spatial distribution. The response of a spectrograph can 
be built in quite the same way, apart the explicit wavelength 
dependence of the monochromatic PSF location. 
PSF construction is therefore a quite general tool. 

An accurate PSF description for an optical instrument, 
either in design or observation phase (through calibration), 
is crucial to any astrophysical measurement. In the former 
case, system parameters are tuned to improve the perfor- 
mance, whereas in the latter the instrument signature must 
be identified for measurements correction and determination 
of the intrinsic objects characteristics. 

The PSF is usually computed numerically, although useful 
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analytical approximations are derived e.g. in iBraati l|2002l) . 
The FFT is used for PSF construction also by professional 
ray tracing packages (e.g. Zemax or Code V). 
We discuss the numerical implementation of the FT, through 
a method of direct computation which, under appropriate 
conditions, may be more convenient than the usual Fast 
Fourier Transform (FFT) tools. The computing time re- 
quired by the proposed PSF construction method is signifi- 
cantly lower than the equivalent FFT case; also, the result 
precision can be higher, since the generation of some arte- 
facts is suppressed. Both aspects are relevant to applications, 
with benefits on processing, resolution and/or precision. 
In most of our simula tions, the case study i s the telescope 
of the Gaia mission (|Perrvman et al]|200ll ). approved by 
the European Space Agency for launch in 2011, and aimed 
at a high precision astrometric survey of our Galaxy. The 
science data simulations require high accuracy, for proper 
validation of the reduction algorithms under development 
by the astrophysical community involved. A specific ap- 
plication to wavelength dependence, a lready dealt with by 
th e authors (|Gai and Cancellierelliooa ). is further developed 
in ICancelliere and Gail ( 2007 ). recently submitted. We also 
mention the possible application of the proposed PSF con- 
struction method to the case of Extremely Large Telescopes, 
which among other technological challenges also feature im- 
pressive numerical requirements. 

In section [2] we recall the basics of the diffraction integral 
and discuss its key implementation features. In section [3] we 
compare some of the more relevant features of FFT and DFT 
implementations, in terms of processing time and accuracy. 
Then, in section f?] we discuss the PSF dependence from 
wavelength and the possible introduction of artefacts. In sec- 
tion [5] we assess the application of the proposed method to 
modern Extremely Large Telescopes. Finally, we draw our 
conclusions. 
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2 THE DIFFRACTION INTEGRAL 

The definitions related to imaging performance derive from 
a wide range o f applications; hereaft er, we mostly follow the 
notation from lBorn and WolJ l| 19851 ). 

A point-like source at infinity, as can be considered a star 
except for very high spatial resolution instruments (like in- 
terferometers), generates a spherical wavefront which can 
be considered flat at the entrance of a telescope. This de- 
scription neglects the disturbances in the medium, e.g. air 
turbulence; in modern ground based telescopes, an Adaptive 
Optics (AO) system is used to recover, at least in part, the 
ideal diffraction limited performance. 

The telescope, ideally, folds the input wavefront towards 
a single point of the focal plane (FP), where the diffraction 
image is achieved. The imaging performance of a real optical 
system can usually be described by the discrepancy with 
respect to ideal propagation, mathematically represented by 
an equivalent deviation from the input flat wavefront, i.e. the 
wavefront error (WFE). 

The diffraction image from a real telescope is built by the 
square modulus of the diffraction integral, here represented 
in circular coordinates {r, <j>} and {p, 9}, respectively on 
FP and pupil: 



I(r,<f>) = k 



pdp d6P(p,9)e 



-iirrp cos(0 



(i) 



The integration domain corresponds to the pupil; for a cir- 
cular aperture, 0^/9^1; 0^0^ 2n, where p is the 
normalised radius. The constant k provides the appropri- 
ate photometric result associated to the source emission, 
exposure time and collecting area; other scaling factors are 



neglected. The pupil function P (p, 6) 



depends 



on the phase aberration function $, which is usually ex- 
panded in a series of terms, e.g. the Zernike functions (j} n 
|Born and Wolj|l985h : 



$(p,(9) = —WFE- 
A 



2tt \ - 



A n (j) n {p,6) 



(2) 



The WFE itself is independent from wavelength, contrarily 
to the pupil function, which is affected by the 2ir/\ factor. 
Also, the nonlinear relationship between the set of aberra- 
tion coefficients A n and the image is put in evidence by 
replacement of Eq. Q in Eqs. [Tlor lA2l 

The non aberrated case, corresponding to WFE — 0, 
can be solved, providing an analytic representation of the 
diffraction image of simple systems recalled in appendix lAl 
The form of the diffraction integral, used in Eqs. [T] and IA2I 
leads naturally to its implementation based on FT tech- 
niques. 



2.1 Discrete and Fast FT 



The D iscrete Fourier Transform (DFT), e.g. from iBracewelll 
(2003), of a one-dimensional sequence s n = s(z n ) of N sam- 
ples of the signal s, function of a time variable z, is: 



JV-l 
n=0 



(3) 



where uj is the pulsation, conjugated to the signal argument. 
The signal s is supposed to be smooth and with finite du- 



ration, with N contiguous non-zero samples; this ensures 
proper convergence of Eq. [3] S(ui) is evaluated on a set of 
M points in the desired interval. Apart variable transfor- 
mation, the diffraction integral in Eq. [T] actually appears as 
the FT of the pupil function, with the natural extension to 
two dimensions of Eq. [3] and the discrete numerical form 
required for practical computation cases. Some remarks on 
sampling resolution and processing time are provided inlBl 

In the Fast FT (FFT) algorithm, significant improve- 
ments in the computation time are achieved under the re- 
striction (usually quite acceptable) that the signal is uni- 
formly sampled, so that s n = s(nT), with sampling period 
T. The computation of all N transform values, for each u) q , 
is performed as a single process, split in a hierarchical se- 
quence of lower order FFT steps. 

The diffraction integral is physically limited to the real 
pupil size; however, with this limitation, only two points 
are generated over the characteristic length of the system, 
which is clearly not sufficient to provide an acceptable image 
detail. A typical solution consists in a formal extension of the 
pupil domain to larger size, allowing to achieve the desired 
resolution, provided the pupil function P is set to zero in the 
external region. In practice, the real pupil must be extended 
to a dummy pupil M times larger, in order to force the 
computation of M points in the Airy radius. The idea of 
zero padding for specific FT computations h as been applied 
e.g. to the packing theorem IjBracew ell 2003); it is a common 
approach also in ray tracing packages (Zemax, Code V). 

The pupil resolution must be sufficient to sample the 
intrinsic variation of the WFE, e.g. placing Q points over 
the pupil size. This generates an array of N = M x Q points 
mapping the dummy pupil along one dimension; notably, 
(M—1)Q of them have zero value. This array can be mapped 
by DFT onto M x Q points over the FP, and the transform 
array covers Q times the Airy radius. 



2.2 PSF construction by FT 

In the PSF construction according to the above prescription, 
the FFT approach involves a significant amount of compu- 
tation over quantities on pupil set to zero. Also, the FP 
region of interest (ROI) is usually not as large as implic- 
itly computed by the FFT, because after a few Airy di- 
ameters the PSF becomes vanishingly small for most op- 
tical systems. This consideration led us to the alternative 
approach of direct DFT computation according to Eq. [3] 
l|Cancelliere and Gail 120071 ) . The key features of this strat- 
egy are: 

- restriction to the physical pupil only; 

- restriction to a limited FP region. 

Thus, the case considered concerns the comparison between 
FFT applied to a square image of logical format N x N, 
where N — M x Q, and DFT applied to an input array 
Q x Q to generate an array M x M. For convenience of some 
numerical examples below, we will set M = Q, although this 
is by no means a restriction to the method. 



3 TEST OF FFT AND DFT 

The first relevant comparison between FFT and DFT al- 
gorithms concerns their effectiveness, i.e. their capability 
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Figure 1. RMS discrepancy between FFT, DFT and analytic 
non aberrated PSF vs. number of sampling points. 
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Figure 2. Mean discrepancy between FFT, DFT and analytic 
non aberrated PSF vs. number of sampling points. 



of providing the correct result. Both are derived from the 
same diffraction integral, so that, apart the case of ac- 
tual coding errors, they both rely on the same conver- 
gence conditions, and are expected to provide the same 
result, within numerical errors. The simulation is imple- 
mented on a desktop computer with a Pentium II proces- 
sor, cl ock 2 GHz, 512 MB RAM , in the Matlab environ- 
ment (http://www.mathworks.com). The Matlab fft func- 
tion is based on the FFTW libr ary (h ttp://www.fftw.org), 
described in lFrigo and Johnson! l| 19981 ). which is commonly 
considered to be quite efficient, whereas the DFT is imple- 
mented by straightforward coding of Eq.[T]for the rectangu- 
lar aperture case, as in appendix |A1 



3.1 Non aberrated PSF 

The first test is performed on non aberrated images, since 
this is the only practical case in which an analytic represen- 
tation of the diffraction image is readily available (appendix 
[A")) . The Gaia telescope parameters are used: rectangular 
aperture, size D x = 0.5 m times D y — 1.4 m (high resolu- 
tion in the y direction); focal length F = 35 m. The ROI on 
the FP is restricted to ±600 u.m (x) x ± 200 (im (y), with 
resolution respectively 6 and 2 u.m, and the DFT is com- 
puted only over this region, using the FFT sampling points 
over both pupil and FP. 

The DFT or FFT precision improves with increasing 
resolution of the sampled function, i.e. with increasing num- 
ber of points. The PSF discrepancy between the FFT and 
DFT result, respectively, and the analytic expression, is 
computed vs. increasing resolution in the pupil plane, over 
the range Q = 10 to 64 sampling points. Since the ac- 
tual FFT image format is Q 4 , the largest array generated is 
iK x 4K, i.e. 16 million pixels. At this point, our desktop 
computers already have a significant virtual memory usage. 
The discrepancy is evaluated as mean, peak to valley (here- 
after "peak for simplicity) and RMS value of the difference 
between computed PSFs, evaluated on the whole ROI. The 
RMS discrepancy of both FFT and DFT with respect to 
the analytic PSF (shown in Fig. [1] respectively by dots and 
circles) decreases with increasing sampling resolution, as ex- 



pected; also, the RMS difference between FFT and DFT 
(dashed line) is quite small, ~ le — 15 of the PSF peak. 
Since the minimum real number used (the Matlab constant 
"eps") is 2.22e- 16, the discrepancy between DFT and FFT 
results is quite close to the limiting numerical noise. Thus, 
both FFT and DFT compute the desired approximation to 
the analytic PSF, with negligible difference between them. 
A similar improvement with sampling resolution is shown 
by the peak discrepancy of FFT and DFT vs. analytic PSF; 
also, the mean discrepancy remains quite small (order of 
le — 15 of the PSF peak value) over the whole range, as 
shown in Fig. In particular, when the pupil is sampled 
with more than 14 x 14 points the RMS discrepancy over 
the ROI drops below 0.1% of the peak value, and the le — 4 
level is reached for more than 45x45 points in the pupil. Sim- 
ilarly, the peak discrepancy is below 1% when using more 
than 14 x 14 points and below 0.1% with sampling above 
48 x 48 points. 

3.2 Aberrated PSF 

The non aberrated PSF case, although quite useful for com- 
parison of the two FT methods, is not particularly interest- 
ing from an optical standpoint, since any real system has 
deviations from the ideal case. However, for a general aber- 
rated case, an analytic representation is not available, so 
that the above simulation cannot be extended to the images 
related to an arbitrary case. However, we can build a set of 
aberrated cases and evaluate the PSF discrepancy between 
DFT and FFT results, which can be compared with the cor- 
responding non aberrated case. 

We generate a set of 1000 aberration cases, with sampling 
resolution of 24 points, using a Gaussian distribution of the 
Zernike coefficients with a — 100 nm. The resulting RMS 
WFE distribution is shown in Fig. [3] The average value is 
94 nm, i.e. an acceptable value for realistic optics in the vis- 
ible (~ A/6 at A = 600 nm), with standard deviation 36 nm. 

For each set of aberration coefficients, the PSF is built 
by both FFT and DFT, then the FFT result is cut to the 
size of the selected FP ROI, where also the DFT is com- 
puted. The discrepancy between the two PSF versions is 
evaluated, deriving for each case the mean, peak and RMS 
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Figure 3. Histogram of the RMS WFE Figure 4. Processing time vs. sampling points of FFT (stars) 

and DFT (crosses), and corresponding fitting functions. 



Table 1. Statistics of mean, peak and RMS discrepancy between 
DFT and FFT result, over the sample of 1000 aberration cases. 
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- 15 


1.9e - 


- 15 


ROI RMS discrepancy 


1.5e- 


- 16 


l.le - 


- 16 



values over the ROI, and recording the values for the whole 
the statistical sample of aberration cases. 

The results are shown in table [U The PSF difference 
between the FFT and DFT methods, for a general case of 
aberration, remains quite small, comparable to their discrep- 
ancy for the non aberrated case (~ le — 15 of the peak), and 
compatible with the numerical noise already experienced in 
the previous test. Therefore, also for the general case of an 
arbitrary PSF, the DFT method generates a very good rep- 
resentation of the FFT result. 

3.3 Processing time 

The processing time for a one-dimensional FFT on N points 
is order of NlogN, and ~ N 2 for DFT. We deal with two- 
dimensional images of logical format N x N = Q 2 x Q 2 , in 
which the region of interest is of order of Q x Q pixels. The 
FT computational load can be explicitly derived for the bi- 
dimensional case corresponding to PSF computation. 
We consider the case in which the one-dimensional FT al- 
gorithm is repeated N — Q 2 times for the FFT, and Q 
times for the DFT, although this is not necessarily the most 
efficient implementation. Then, the cost for bi-dimensional 
FFT and DFT becomes of order of TV 2 log iV 2 = Q 4 log Q A oc 
Q 4 logQ, and Q 2 x Q = Q 3 , respectively. Introducing ad- 
ditional numerical factors (e.g. the Airy diameters in the 
ROI), the scaling law with respect to the number of sam- 
pling points Q does not change, apart the coefficient. 
The potential gain of DFT vs. FFT is thus ~ QlogQ. 

Notably, similar gain may be expected in other ap- 
plications in which the input data or the FT results are 
required with high resolution only over a restricted re- 
gion for both conjugate variables. We remark also that the 



condition of uniform sampling is not necessarily true for 
any measurement, and can be sometimes inconvenient. The 
above estimate of processing time is valid only within the 
computation-limited regime, and for any real computer, at 
increasing size of the processed arrays, the processing be- 
comes input/output-limited, when the physical memory is 
saturated and the virtual memory mechanisms start swap- 
ping data towards the mass storage devices. 

A stopwatch timer function of Matlab is used to keep 
track of the execution time of both FFT and DFT. The sim- 
ulation is repeated for 10 iterations, in order to reduce the 
fluctuations due to operating system overheads. The pro- 
cessing time is averaged over the 10 sets, providing the re- 
sults shown in Fig. [4] where the fit to the analytic functions 
<2 4 logQ (FFT) and Q 3 (DFT), is also shown. The fit to 
experimental data is reasonable, although noisy; the fit pa- 
rameters are respectively 129 ns (FFT) and 548 ns (DFT). 

Over the range of cases considered, DFT is significantly 
more efficient than FFT, and the growth rate of the compu- 
tational cost is lower. 



4 WAVELENGTH DEPENDENCE 

Since in our examples the FP resolution is fixed (2 u.m), the 
pupil resolution scales with the wavelength A: the pupil is 
better sampled at shorter wavelength and vice versa. The 
variable pupil resolution is potentially associated to arte- 
facts in the images, depending on wavelength, leading t o 
potential measurement errors, as in lBusonero et all |2006). 
In order to test this aspect, the wavelength range Ai = 
590 nm to A2 = 610 nm is explored with AA = 1 nm resolu- 
tion. The case considered is the non-aberrated optical PSF, 
because of the known analytic representation. We evaluate 
the discrepancy with respect to the analytical function of 
three different PSF representations, respectively built one 
by FFT, and two by DFT, the former using the FFT pupil 
sampling resolution (DFT 1), and the latter using fixed sam- 
pling of 64 points on the pupil (DFT 2). Notably, DFT 2 
resolution is lower than FFT or DFT 1, using respectively 
136 (77) and 146 (£) points at A = 600 nm. 

The variation with wavelength of the three cases of PSF 
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Figure 5. PSF discrepancy vs. wavelength of FFT (circles), 
DFT 1 (dots) and DFT 2 (stars). 
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discrepancy with respect to the analytic representation is 
shown in Fig. [5] respectively for the results from FFT (cir- 
cles), DFT 1 (dots), and DFT 2 (stars). The DFT 1 resolu- 
tion was chosen to match the FFT resolution, so that it is 
reasonable to expect similar behaviour. We remark the large 
variation in discrepancy of FFT and DFT 1 over a compara- 
bly narrow spectral range. This is expected to lead to poten- 
tially significant errors in the construction of polychromatic 
images, depending on the selected spectral sampling. 
Besides, using a constant pupil sampling vs. wavelength (as 
for DFT 2, compatible with fixed FP resolution), we get a 
uniform distribution of PSF discrepancy vs. wavelength. In 
spite of lower pupil resolution, the precision of DFT 2 is 
significantly better over most of the wavelength range, and 
above all it is constant, i.e. wavelength independent. 



5 EXTREMELY LARGE TELESCOPES 

The next generation of telescopes currently studied aim at 
a significant improvement in sensitivity with respect to cur- 
rent instrument ation, as for the O verwhelmingly Large tele- 
scope (OWL), in lGilmozzil (|2004l ). Of course, the related cost 
and engineering challenge is also significant. Other proposed 
Extremely Large Telescopes (ELTs) also have typical size of 
a few ten meters, sometimes achieved as diluted apertures, 
as proposed e.g. in lMartin et alj l|2003l ). for a better trade-off 
between complexity and sensitivity. 

The PSF computation for such instruments is in prin- 
ciple based on the same tools (diffraction integral) used for 
previous telescopes. Besides, processing improvements may 
be appealing to the purpose of increased precision in the 
analysis and more efficient design development. 

In recent propositions, the OWL diameter is 42 m; as 
an example, we assume operation in the K band (2.2 urn), 
corresponding to an Airy diameter of 26 milli-arcsec, and 
requiring an effective focal length of 560 m for its imag- 
ing on four 18 u.m pixels. For simulation of the OWL PSF, 
we assume a minimum sampling of order of 10 points per 
pixel, i.e. 40 points over the Airy diameter, corresponding to 
1.8 urn resolution on the focal plane. We set again the ROI 
size on the FP to 6 Airy diameters (on each dimension), for 
a logical format of 240 x 240 pixels. 



Figure 6. Example of an OWL PSF with static aberrations and 
AO residuals. 

For proper sampling of the atmospheric turbulence over 
the telescope pupil, in K band, we set a resolution of order 
of 0.1 m, requiring an array of 420 x 420 points on pupil. 
However, the FFT size associated to the problem, due to the 
compound effect of 40 points in the Airy diameter times 420 
points on the pupil, is 16800 x 16800, simply unmanageable 
on the desktop class computer used in our simulation. 

The simulation is set with a few Zernike coefficients set 
to non zero values, to mimic a realistic optical response, and 
with a Gaussian random residual from AO correction, with 
a — 300 nm. The corresponding RMS WFE on the pupil in 
the case shown is 310 nm. The resulting PSF is shown in 
Fig. [5] in logarithmic units; the signal profile is affected by 
perturbations both on low and high frequency, respectively 
from optical aberrations and AO residuals. The monochro- 
matic PSF computation requires about 0.6 s, increasing to 
about 5 s for 0.05 m pupil resolution. Therefore, analyses 
of such optical system, either with respect to instrumental 
aspects (including AO response) or to astrophysical perfor- 
mance, can be done without high-end computers. Besides, 
for a given computing power, the proposed method allows 
a significant increase of the range of study cases, e.g. for 
statistical evaluation. 



6 CONCLUSION 

When PSF computation is required only on a limited region 
of the focal plane, DFT computation can be significantly 
more efficient than FFT, providing the same results (within 
the numerical noise) with much lower processing time. This 
is because computation is restricted to the actual points re- 
quired by the problem, rather than those imposed by the 
FFT definition on both input and output variables. Both 
FFT and DFT produce a good PSF representation, with pre- 
cision increasing with the number of sampling points used. 
Artefacts can be introduced in the PSF, e.g. as a function of 
wavelength, depending on the sampling strategy; DFT pro- 
vides more flexibility and can thus improve the precision. 
Processing time is relevant also with respect to preci- 
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sion of the PSF computation and of the performance anal- 
ysis based on such data; in particular, precision is in many 
cases dependent on the number of sampling points used for 
representation of the relevant problem variables. Therefore, 
the capability of processing higher resolution data, in a given 
time frame, ensures higher precision on each PSF; besides, it 
is also possible to take advantage of the proposed method to 
improve the statistical sampling of the problem domain, by 
generation of a larger amount of simulated data, depending 
on the most convenient trade-off. 

Also, the development of new generation, large size tele- 
scopes may benefit from faster and more accurate computa- 
tion tools, for both design and science assessment. 
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APPENDIX A: 
IMAGES 



IDEAL DIFFRACTION 



In the simple case of an unobstructed circular pupil of di- 
ameter D, at wavelength A, and with focal length F, we get 
the radial symmetric PSF of Eq. IA1I 



I(r) 



Jx (jrrD/XF) 
TirD/XF 



(Al) 



Here Ji is the Bessel function of the first kind, order one, 
and r the radial FP coordinate. The image characteristic 
size is the Airy diameter: 2A4XF/D. 

For a rectangular pupil, it is convenient to use cartesian 
coordinates, e.g. {x, y} on FP and {£, r)} on the pupil, with 
integration performed over the appropriate region: 
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Figure Al. Non aberrated PSF of Gaia: logarithmic intensity 
at 600 nm. 



I(x, y) = k 



(A2) 



As an example, we consider the current Gaia telescope ge- 
ometry, i.e. a rectangular aperture, with size D x = 0.5 m x 
D y = 1.4 m (high resolution in the y direction) and focal 
length F = 35 m. The PSF is described by the sine squa red 
function, where sine z = (suittz) / (tvz) (|Bracewelfll2003h : 



( xDx \ 
\XF ) 



( yDy 

\ \F 



(A3) 



A representation of the non aberrated PSF, at A = 600 nm, 
is shown in Fig. IA1I in logarithmic units normalised to the 
central peak intensity. 



APPENDIX B: 
PROCESSING 



FT RESOLUTION AND 



Usually, a uniform sampling period T is assumed for DFT, 
although this is not strictly implicit in the definition [3] The 
signal transform is computed at pulsation u>„ = 2nn/ (NT), 
where n = . . . (N — 1), so the resolution is Alj — 2n /(NT) 
(or Av = 1/(NT), in terms of frequency: u) — 2nu). The 
DFT thus associates the resolution in one domain to the 
full range of the corresponding conjugate variable. 

The computational load of FFT and DFT of a sequence 
of TV values is respectively of order of TV log TV and TV 2 , so 
that for increasing TV the FFT is much more convenient of 
DFT performing the same computation. The expression is 
valid only when the number of samples is a power of two: 
N — 2 n , for some integer n; for other values, the result is 
highly dependent on the selected algorithm. 

The numerical implementation of the DFT for the rect- 
angular pupil case follows from the definition, apart de- 
tails of proper definition of the pupil and FP coordinates, 
respectively {x n , y m } and {£ p , r/q}, with suitable indices 
{n, m, p, q}: 
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